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, We describe an implicit procedure for solving linear equation systems resulting from the discretiza- 

tion of the three dimensional (seven variables) linear Fokker-Planck equation. The discretization 
of the Fokker-Planck equation is performed using a twenty-five point molecule that leads to a co- 
efficient matrix with equal number of diagonals. The method is an extension of Stone's implicit 
procedure, includes a vast class of collision terms and can be applied to stationary or non station- 
ary problems with different discretizations in time. Test calculations and comparisons with other 
V-^ , methods are presented in two stationary examples, including an astrophysical application for the 

Miyamoto-Nagai disk potential for a typical galaxy. 

r- ; 

SpH 1 I. INTRODUCTION 

In dealing with solutions of partial differential equations we often encounter a set of linear equations that has to 
Q ■ be solved. This set of linear equations depends on the method used for discretization. In general, when dealing with 
three dimensional systems the number of linear equations increases and the numerical solution of these equations 
uses most of the computing time. An extreme case is the Fokker-Planck equation. The Fokker-Planck equation is 
also known as the Fokker-Planck approximation because truncates the BBGKY (N.N. Bogoliubov, M. Born, H.S. 
Green, J.G. Kirkwood, and J. Yvon) hierarchy of kinetic equations at its lowest order by assuming that correlation 
between particles only plays a role as a sequence of uncorrelated two-body encounters 0, Q|. Note that the only 
| "approximation" made in the Fokker-Planck equation comes from the model adopted for collisions and, in fact, the 
i '"T 1 - Fokker-Planck equation can be derived from first principles and no ad hoc suppositions are needed. The solution of 
the Fokker-Planck equation equation is not an easy task because in the three dimensional case it has seven variables: 
three space coordinates (x), three velocity coordinates (v) and time (t). In the two dimensional case there is a 
simplification because the total number of variables are five. In either case, the large number of grid nodes needed for 
the computation of the solution becomes a storage data problem. In a three dimensional problem, for a stationary or 
non stationary equation, the number of linear equations corresponds to the number of nodes in the phase-space grid 
iy-j \ ( X 5 V )- If we divide each of the phase-space variables' interval of the distribution function in nine parts (ten nodes), 
we will have a grid with 10 6 nodes. In a simple numerical method we have to store and solve a matrix with 10 12 
' elements. For the two dimensional case, the main matrix will have 10 8 elements. With 10 grid nodes per variable 
only very simple geometries can be described. The large number of matrix elements brings us another computational 
problem, the slowness of the codes. In the discretization process of the Fokker-Planck equation, a system of linear 
equations is obtained and arranged into a matrix form (coefficient matrix). For the case of a finite difference scheme 
discretization in three dimensions with a twenty five point molecule, we see that approximately less than 0.003% of 
the elements are different from zero. This incentives us to search for alternative and faster methods, usually iterative, 
>^ to solved the linear system using only the non null data. Note that, in general, the coefficient matrix is not symmetric. 
Ft" 1 ! So, powerful methods like the Conjugate Gradient 0, Q and Cholesky Q| decomposition can not be used. Our main 
J> . goal is to obtain a code that allows us to obtain a fast and effective numerical solutions on high resolution schemes of 
J>Tj ' the three dimensional linear Fokker-Planck equation in a direct way. The importance and difficulties of having three 
dimensional solutions of the Fokker-Planck equation can be summarized in the words of Binney and Tremaine (p| 
page 245), here in relation with galactic dynamics: Finding the particular function of three variables that describes any 
given galaxy is no simple matter. In fact, this task has proved so daunting that only in the last few years, three-quarters 
of a century after J cans' s (Jeans |||) paper posed the problem, has the serious quest for the distribution function of 
even our own Galaxy got underway. 

We mean by direct numerical calculations of the Fokker-Planck equation a method that is neither statistical nor 
mean field approximation Numerical solutions can be performed using statistical approximate methods like the 
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method of moment equations 0jHi@ an d Monte Carlo methods El EI- Also, solutions have been found using the 
orbit-averaged Fokker-Planck equation with action-angle variables [f| , this last method reduces the equation involving 
six phase-space coordinates plus time to one involving only three actions plus time. Another method that have been 
used for solving the Vlasov equation with good results is the operator splitting technique [T2I Il3l ] . Two dimensional 
numerical integration of the Vlasov equations can be found in [l4j |. The problem associated with this method is 
that, even in the two dimensional case, the time spent in making the interpolation is very expensive. Direct three 
dimensional numerical calculations of the Fokker-Planck equations, to the best of our knowledge, has not been done. 
As we said before, the main difficulty to find the solution of the three dimensional Fokker-Planck equation is the large 
number of linear equations that is translated in a high computational cost involve in the process. In this article we 
present a variations of Stone's method that leads us to solve with low computational cost the three dimensional 
linear Fokker-Planck equation. Variations of Stone's method has been applied to other situations in two ^(| and 
three |l7l | space dimensions when dealing with fluid flow, heat transfer and Laplace equation. 

In this article we study how to solve the system of equations that arises from the discretization process using a 
finite difference scheme, in particular we used a twenty-five point molecule, but we have to keep in mind that other 
discretization procedures gives practically the same form for the coefficient matrix and this method can also be applied, 
i.e. sparse matrices with the same number of diagonals. For example, in two dimensions, the diffusion equation can 
be discretized in an uniform rectangular grid by using the finite difference scheme, the finite volume method (which 
is a discretization of the equation in integral form) or the finite element method, see |T^| . All of these methods gives 
five diagonals, in the finite element method the number of diagonals depends on the type of interpolation considered, 
and the coefficient matrix can be solved by the same numerical method. So, in this article, we shall not discuss what 
discretization method is better for the problem considered, we rather present a numerical method that allows us to 
find the solution of a coefficient matrix with twenty-five diagonals that can be obtained with different discretization 
methods. 

The article is organized as follows. In Section II we present the linear Fokker-Planck equation to be solved. We 
consider a linear collision term that includes a vast class of collisions. In Section III, we describe the algorithm of 
the modified Stone method to obtain the incomplete LU decomposition. The discretization of the Fokker-Planck 
equation is performed using a central difference approximation for the phase-space variables (x,v) that is described 
with a twenty-five point molecule. This point molecule also allows to describe different discretizations in time, as the 
implicit Euler or Crank- Nicolson discretization. The derivation of the LU decomposition is made without assuming 
any particular discretization in time, maintaining its derivation as general as possible. The notation used for the 
diagonals in the coefficient matrix is also shown. In Section IV, we test our algorithm with an astrophysical example 
by solving the Fokker-Planck equation for the widely used Miyamoto-Nagai disk potential 0] in three dimensions. 
The parameters of the model are chosen to represent the Newtonian potential of a typical galaxy. Also, we used 
a Fokker-Planck test equation to compare our algorithm with the Generalized Minimal Residual Method (GMRES) 
that solves large sparse matrices. In Section V, we show how to modified the code to implement curved boundary 
conditions. Finally, in Section VI, we summarized our results. 



II. THE GENERAL PROBLEM 



The linear Fokker-Planck equation can be written as 



^ + v-v/+v-v„/ = r[/], (i) 

where v represent the velocity of the particles, V is the usual gradient, V„ is the velocity gradient (derivations are 
done with respect to the velocities), v is the acceleration and the symbol T[f] denotes the rate of change of / due to 
encounters (collision term). We consider as the collision term the expression 



r[/] = A(x,v)V 2 / + S(x,v)V2/ + C(x,v)V/ + i?(x,v)V u /+ ]T ^(x,v)^-, (2) 

where A, B, C, D and Eij are arbitrary functions of the phase-space variables x and v. The equation above describes 
a vast family of collisions. In particular, with the mixed velocity derivative term present in we can take into 
account the important collision term found by Rosenbluth et al. [2(| used in gravitating systems and plasma physics, 
see for example 0, and [2II. |24| respectively, and reference therein. If we need mixed space derivatives instead 
of mixed velocity derivatives, we can use the same code presented in this article to solve the problem. If a particular 
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TABLE I: Nomenclature and relations between the matrix form and the one dimensional storage index at node p of the 
twenty-five terms used in the discretization of the Fokker-Planck equation. 
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problem requires the inclusion of mixed velocity derivatives as well as mixed space derivatives, it is possible to develop 
a similar numerical procedure following the steps of this article, but it complicates the incomplete LU decomposition 
used for the method, i.e. we need twelve extra diagonals on the coefficient matrix (thirty-seven instead of twenty-five 
diagonals). 

In Section IV we solve numerically the stationary linear Fokker-Planck equation (I}. In general, Eq. is non- 
linear because v = F(/)/m, where F is the force and m is the mass of the particle. Another kind of non-linearity may 
arrise from the collision term considered. A way to deal with this kind of problem is to start with a given distribution 
function at time t from which we can calculate the force and collision term at this time. Then, this force and collision 
term are replaced into the Fokker-Planck equation from which we obtain the distribution function for a later time 
t + At. With the recently calculated distribution function we can calculate again the force and collision term at time 
t + At and the process is repeated. An application for a stationary non-linear problem can be found in [25|. in which 
we found the distribution function that satisfies both Fokker-Planck and Poisson equation in two dimensions for a 
Kuzmin-Toomre thin disk. 



III. DESCRIPTION OF THE ALGORITHM 

The system of equations obtained from the discretization of the Fokker-Planck equation (IIIL'I) using the central finite 
difference approximation for the phase-space and a temporal discretization in time (implicit Euler, Crank-Nicolson, 
etc.) can be cast (for each time step) into the simple form, 



A* = Q, 



(3) 
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FIG. 1: Form of the coefficient matrix A obtained from the discretization of the partial difference equation of our problem. 
The separation between the main diagonal P and the other diagonals are indicated. Also, the distance between the diagonals 
of the mixed derivatives to their nearest diagonals are shown. Note that the figure is not in scaled. 



in which A is a square coefficient matrix N noc ie X N no de (N no d e the number of nodes in the discretization grid), W 
is the vector matrix of the nodal variable values, and Q is the source vector. The position of the grid nodes in the 
phase-space (x, y, z, v x , v y , v z ) is performed by six indexes (i, j, k, I, m, n), where i represents the index for the variable 
x, j represents the index for the variable y, etc. The ordering of nodes in this six dimensional space is made as 
follows. The surface n=constant are stacked one above another. Within the fifth dimensional space (for each n) the 
hyper-surfaces m=constant are stacked one above another. Within the fourth dimensional space (for each n and m) 
the hyper-surfaces Z=constant are stacked one above another. Within the three dimensional space (for each n, m and 
I) the surfaces fc=constant are stacked one above another. Within the two dimensional space (for each n, m, / and 
k) the index j increases first (y-direction) than the index i (x-direction) . The one dimensional storage index p of the 
vector matrix \E' is calculated from the six-dimensional grid indexes k, I, m, n), i.e. 

p = (n - l)n ijklm + (m - l)n, ljk i + (I - l)n ijk + {k- l)n l:j + (i - l)rij + j, (4) 

with 

i = 1 • • • rii ; j — 1 • • • rij ; k = 1 • • • rik ; I = 1 • • • ni ; m = 1 • • • n m ; n = 1 • • • n n ; 

n-ijkim — nirijn k nin m ; n ijk i = riinjn k ni; n i j k =n i njn k \ = niUj] (5) 

where rij, rij, n k , ni, n m , n n denote the number of grid points for each variable. Therefore N no( i e = ninjn k nin m n n . 
With the help of the storage index p we can switch each point of the twentyfive-point molecule from the matrix form to 
the one dimensional position representation. This is done in TableQJby making the equivalence f(x, y, z, v x , v y , v z ) = 



FIG. 2: Schematic representation of the matrices L, U and their product M. The multiplication of the matrices L and U leads 
to extra diagonals (dotted lines) not present in the coefficient matrix A. The diagonals of L, from the left-bottom corner to the 
main diagonal, are: DD3, DB3, D3, DT3, DU3, DB2, D2, DT2, Dl, W, S, P. The diagonals of U, from the main diagonal 
to the right-up diagonal, are: 1, N, E, T, Ul, UB2, U2, UT2, UD3, UB3, U3, UT3, UU3. The diagonals in M are products 
of these two sets of diagonals but we have to be careful because more than one product can be at the same diagonal. 

^>(i,j, k, I, m, n). Until now, we considered only the non stationary case, but the discretization in Table[I]can be used 
in the non stationary as well as the stationary Fokker-Planck equation. In both cases, the final system of equations 
is of the form @, with A being a sparse matrix with elements different from zero in only twenty-five diagonals, see 
Fig. ^ Now, we want to develop an iteration method in order to solve the system of equations Q . After h iterations 
of such method, the approximate solution ^ h do not satisfies © exactly, their is a non zero residual R such as 

A^ h = Q-R h . (6) 

The purpose of the iteration procedure is to drive the residual term to zero after some number of iteration (actually we 
stop the iteration when the residual term attained some imposed small value condition). Let us consider an iterative 
scheme for a linear system in the form 

M^ h+1 = 0^ h + B, (7) 

when convergence is achieved we must have that A = M — O and B = Q. An alternative version of this procedure 
can be obtained by subtracting M^ h from both sides of to have, 

A/A' l+1 = R h , (8) 

where A' l+1 = - and R h = B - (M - O)^' 1 = Q - A$ h . Any effective iterative method to solve must be 
cheap and converge rapidly. For faster convergence, the matrix M have to be a good approximation of the coefficient 
matrix A, i.e. we must have 0^! h small. The original idea of Stone is to use for the iteration matrix M an incomplete 
LU decomposition of the matrix A. The reason for this choice is that LU decomposition is an excellent linear system 
solver. The matrices L and U have elements different from zero only in the diagonals in which A have also elements 
different from zero. The product of the matrices, L and [/, provide a matrix M with a larger number of diagonals 
with elements different from zero, see for instance Fig. |5J To make the decomposition LU unique, we set the elements 
on the principal diagonal of U equal to 1. In doing the LU multiplication, we have to pay extra attention because 
sometimes more than one diagonal product are at the same distance from the main diagonal P, e.g. the product 
Ld3 • U± and Ldd^ ■ Ujji are at the same diagonal in M. The multiplication rules for matrices furnish the elements of 
M = LU at node p (see Appendix |A"|> . where for example, Mw\uu3 represents the diagonal in M that is obtain from 
the multiplication between the elements of the diagonal S in L with the elements of the diagonal UU3 in U. 

Now, we choose L and U in such a way that M(= A + O) is the best possible approximation to A. The standard 
method for decomposition is to let O to have elements different from zero in the diagonals of M that corresponds to 
the diagonals not present in A and to force the other diagonals of M to be equal to the corresponding diagonals in A. 
But this method converges slowly because the resulting matrix O is not small. Stone recognized that the convergence 
of the method could be faster if we allow O to have elements different from zero also in the diagonals present in A. 
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The key idea is that the contribution of of the diagonals not present in A partially canceled the contribution of 
O^S of the diagonals present in A, in such a way that 

O* w 0. (9) 

Note that in (|Alfl they are diagonals in M that present more than one term. In general, the principal diagonals 
have more than one element, as for example Mp,T3, but beside these principal diagonals there are other no-principal 
diagonals that have more than one element, like M D3 \u B2 - Now, relation © can be written for one grid node in 
several ways. The usual way is to consider the elements of these no-principle diagonals as part of the same diagonal. 
Other way is to consider these elements as they were from different diagonals, thus in this case, following the above 
example, the no-principal diagonal M]J3\UB2 is split into two diagonals Mp,3\uB2 an d Mub3\U2- We obtained the final 
relations for the LU decomposition in both ways and we found that the LU decomposition considering the elements 
as they were from different diagonals is faster by a factor of two. Hereafter we consider this case. The explicit form 
of equation is given in Appendix iBl 

The problem now is to defined the elements of O to satisfy Eq. l|Blfl without introducing additional unknowns. 
If we expect the solution of the partial differential equation to be smooth, we can approximate the values of ^ b\N: 
^n\W: etc, in terms of the values of \& at nodes corresponding to the diagonals of A. Stone proposed the following 
approximations (other approximations are possible), 

^b\n ~ a(^ B + *w -*p), 

^w\n ~ a{^> w + - etc (10) 

where a is a constant. Stability analysis made by Stone requires that a must be between < a < 1. Replacing the 
above approximations into (|B1|I we obtain the elements of O as a linear combination of the elements of M, see for 
instance Appendix[U] Now, using the relation M — A + O together with expressions (|A1|) and (|C1|) . we find that the 
elements of the matrices L and U are given by, 

if = — ^1 [x] = DD3, DB3, D3, DT3, DU3, DB2, D2, DT2, Dl, B, W, S, 
1 + aK[ X ] 

TP _ AP _ tPttP-1 _ TP TjP- n i _ TP TjP- n i] _ TP T rP~ n ijk _ TP T tP~ ("yfel -«i jk ) _ jP TjP-riijkl 
Ijp — lip ^ s N Ij W u E lj B u T ±J D1 U UI ^DT2 U UB2 ^ D2 U U2 

_tP T rP~( n ijkl+nijk) _ jp Tj-p—(nijUm—nijkl) tP jjP—{nijklm—riij k ) _ T p jjP~riij k i m 
lj DB2 U UT2 ^DU3 U UD3 lj DT3 U UB3 ^ D3 U U3 

-^ ra C y ''"" +n " Jt) - LnnM- u i 3 nm,m+nmi) + a(K N + K E + K T + K m + K UB2 + K U2 + K UT2 
+K um + K UB 3 + Kjj3 + K UT 3 + Kuuz), 
Afv-i — 0iK\ x \ — C\ X ] 

UF X] = — , [X] = N, E, T, Ul, UB2, U2, UT2, UD3, UB3, U3, UT3, UU3, (11) 

where the explicit form of the functions Ciji and K\x] ar e given in the Appendix [D] The elements of the LU 
decomposition have to be calculated in the order specified in <|11[) . In doing this, we must take into account that a 
certain element is considered equal to zero if its storage index is less or equal zero, e.g. if p = 3 and rij = 5 then the 
elements with index p and p — 1 are different from zero, and the elements with index p — rij, p — riij, etc are equal 
to zero. When mixed derivatives are not present we must set all the elements with index DD3, DB3, DT3, DU3, 
DB2, DT2, UB2, UT2, UD3, UB3, UT3, UU3 equal to zero. Once obtained the LU decomposition, the system of 
equation is solved combining M = LU with (JSJ) to obtain 

LUA h+1 = R h , (12) 

and here we set 

LT h = R h^ 

UA h+i = T h^ ( 13 ) 

from which we obtain the solution of our problem by solving two triangular systems. In this iterative method, the 
matrix elements of L and U are calculated only once before the first iteration. In other iterations, only the residual 
i?, T and A are calculated using the two triangular system mentioned above, i.e., 
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TP — (JfP — TP TP~( n i3kim+n ijk i) _ TP TP-( n ijki m +nijk) _ TP TP~ n ijkim _ TP TP-( n ijki m -riijk) 
1 — yjx ^UD3 L n DB3 1 n D3 L DT3 1 

TP yp-(n ijk i m -n ijk ,) jP TP-(n ijkl +n ijk ) jP yp-n ijkl jP yp-(n ijkl -n i]k ) _ jP yp-n ijk 
n DU3 1 DB2 1 ^D2 L ^ DT2 1 n Dl 1 

—L p B T p ~ nij - L p w r p - n i - ^T p_1 )/L^, 

^P _ TP _ TjP^^^P+inijMm+riijki) _ jjP^^P+inijklm+iiijh) _ TjP^^P+nijklm _ jyg. ^2>+(«yw m — tty*) 

— TJP AP+( n ijkl m ~n ijkl ) _ jjP <\p+(n ijk i+n ijk ) _ jjP Ap+n ijk , _ tjP \p+(n ijkl -n ijk ) _ jjP \p+n ijk 
U UD3^ U UT2 LA U U2 L - X U UB2 LA U U1 LA 

-U p A p+n » - U p E A p+n > - [7^A p+1 , (14) 



where 



RP = QP — A v p ^p - A p N yP +1 - A^p- 1 - A p E ^P +n > - A p w ^>P- n i - A p T ^ p+n '^ - A p B f p - r ^ - A p ui ^ p+n ^ k 
— A p Dl $> p ~ ni3k — A P JB2 '$iP + ( nijkl ~ nijk ^ — A p DT2 ^ p ~ l > nijkl ~ nijk ^ — A P J2 ^ p+nijkl — A p D2 ^ p ~ niihl 

AP vT/P+( n «fc!+ n yfe) AP \f,p-(nij k i+n ijk ) AP \f,P+(ni jk i m -n ijkl ) AP rt,p-(n ijklm -n iikl ) 
/1 (7T2* /1 DS2* A UC3* A DU3 W 

-A p um ^ p+ ^ kl ™- n ^ - A p DT ^ p - {n ^ klm ~ n ^ k) - A p r3 ^ p+ni]klm - A p ra r"" H ™ 

AP iY/P+(mjklm,+nuk) AP iT/P— (««Wm+n«fc) AP iT/P+(ri«ifclm+n«,-fci ) AP \t,p-(nij k i m +nij k i) 
A UT3^ DB3 ^(7(73* A DD3 * ' 

in which, for simplicity, we have omitted the iterative index h. 



IV. TEST RESULTS AND COMPARISONS WITH OTHER METHODS 



As was mentioned in the Introduction, the huge number of nodes needed to solve the Fokker-Planck equation leads 
to large amount of data that has to be stored in a matrix, this fact reduces the possible codes for testing the results. 
Furthermore, in general, the Fokker-Planck equation is not symmetric, and for that reason other methods can 
not be used. In this case, the Generalized Minimal Residual method (GMRES) l2al is the most appropriate choice. 
The GMRES method belongs to the class of Krylov based iterative methods |27l I28L l2^| and was proposed in order 
to solve large, sparse and non Hermitian linear systems. This section is divided into two parts. In the first part, we 
apply the method presented in this article to a physical example. In particular, we choose the astrophysical problem 
of finding the distribution function of the Miyamoto-Nagai disk in three dimensions. In the second part, we used a 
simplified version of the physical example to make numerical comparisons. We compare our method with the best 
available method in solving huge sparse matrices, i.e. GMRES. We prefer to make the numerical comparisons with 
the simplified version of the physical example because the large amount of parameters present in the physical example 
make the analysis cumbersome. This analysis is important because we show some advantages and limitations of our 
method. 



A. Physical Example 

To test of our code we begin with an important astrophysical application, i.e. solve the Fokker-Planck equation to 
find the distribution function of the Miyamoto-Nagai disk [lJ| in three dimensions in a stationary regime. The Fokker 
Planck equation to be solved is [j| 



v-V/ + v-V l ,/ = -^— [f(x,v)D(A Vi )} + - £ ——lf { ^ y)D{ A Vl A Vj )l (15) 

i— 1 1 ij — 1 % ^ 

where v = — V<I>, $ = GM . = . G is the gravitational constant, M is the total mass of the system, a 

and b are parameters that depending on the choice the potential can represent anything from an infinitesimal thin 
disk to a spherical system, and the functions D(Avi) and D(AviAvj) are known as the diffusion coefficients. These 
diffusion coefficients were calculated by Roscnbluth et al. 20] considering a test star of mass m moving through an 
infinite homogeneous sea of field stars of mass m a who has mean velocity equal to zero. Moreover, the interaction 
between the particles are ruled by an inverse square force, and also, each stellar encounter involve only a single pair 
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of stars and are independent of all others. These diffusion coefficients are simplified if the field stars distribution 
function is a Maxwellian distribution. The explicit form of these coefficients are |5J. 



D(A Vl ) = -D(A V]l ), 
v 



D(Av i Av j ) 



<D(Av 2 ) 



where D(Avn), D(Av?,) and Avj_ are given by 



l -D{Avl)] + ld l3 D(Avl), 



(16) 
(17) 



D(Av\\ 
D(Avf { ) 
D(Avl) = 



^G 2 p{m + m a ) lnAG(X) 



> _ 4\/27rG 2 pm a lnAG(X) 
a X ' 

4V2irG 2 pm a liiA \eri(X) - G(X) 



a 



X 



(18) 
(19) 
(20) 



where p and a are the density and velocity dispersion of the field stars respectively, X = vj (y2a), erf (X) is the error 



function, G(X) 



l 

2 A 1 



erf(X) 



2X „-X 2 



A : 



is the maximum possible impact parameter (usually 



G{m+m a ) ' L 

set of order the radius of the system), and i>t yp is a typical velocity of stars in the system. We shall find the solution of 
Eq. (|15fl in a six dimensional 'box' using ten grid nodes for variable that leads to a problem involving 10 6 unknowns. 
If we want to solve this set of equations using a conventional solver we need to storage 10 12 elements in the coefficient 
matrix. We considered that the ranges of the velocities v x , v y , and v z are between [-350 km s~ 1 , 350 km s _ 1 ] , the ranges 
of the coordinates x and y are between [-40 kpc, 40 kpc], and the coordinate z between [-1 kpc, 1 kpc]. At the borders 



of the 'box' we used a Dirichlet boundary condition of the form / = exp[- 



-y 2 + z 2 )/2a 2 ] 



e X p[-(v 2 x +v 2 y + v 2 z )/2a 



We set the parameters a = 4 kpc, b = 1 kpc, m a = m = M , a s — 10 kpc, a v = 100 km s , u typ = 200 km s , 
bmax = 40 kpc and the total mass M = lO 12 Af . The values set for the parameters correspond to typical galaxy 
data like the Milky Way. Also, as a first approximation we spread the total mass uniformly along the disk and set 
the density p ^constant. In Fig. [21 we present two graphs of the distribution function that represent the numerical 
solution of I jl 5(1 . One is the surface distribution function on the plane x — y and the other is the velocity distribution 



function on the plane v y — v 
Maxwellian form as it should. 



they are plotted at different points of the grid. We note that these figures have 



B. Numerical Comparison 

In the next example we shall compare our code based in the LU decomposition 1(11(1 with the GMRES algorithm. 
The problem of this method is that the storage of the orthonormal basis may become prohibitive for some matrices, 
this storage depends on the value of the restarting parameter. The restarting parameter of the GMRES algorithm 
determine the number of the orthonormal vectors for the Krylov subspace that the code stores in order to calculate 
the updated solution and residual at each time step. At each time step the code stores one vector. After a number 
of steps equal to the restarting parameter, the code construct the most recent update and use it as a first guess to 
restart the next set of iterations. The convergence of the method is guaranteed for large numbers of the restarting 
parameter, but this means that more vectors have to be stored and computer memory problems may appear. We can 
use lower values of the restarting parameter but this increases the time spent in finding the solution. In particular, we 
used the gmres routine of Matlab because it can handle sparse matrices. We also used a public GMRES software [3(| , 
this software allows us to choose between different kinds of preconditioned and orthogonalization procedures but its 
drawback is that can not handle sparse matrices. We choose as a test equation a stationary Fokker-Planck equation 
similar to the physical example of the previous section. Note that in the non stationary Fokker-Planck equation we 
can perform a Crank-Nicolson discretization in time, which is also an implicit procedure, and the method presented 
in this article can be applied. The stationary Fokker-Planck equation considered for the test is 



V/ + v- V„/ = -V„/- 



d 2 f 



dvidvj 



E 



o 2 f 



i—i 1 



(21) 
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Distribution Function 

0.00012 r 





FIG. 3: Numerical solutions of the distribution function of Eq. I|15fl . Top: surface distribution for the case with v x ~ —117 km 
s _1 , v y ~ 39 km s~ , v z ~ —272 km s _1 and z ~ 0.55 kpc. Bottom: velocity distribution for the case with x ~ —22.23 kpc, 
y w 13.32 kpc, z « —0.11 kpc and v x « 117 km s~ . 



where v = — V<I>, $ = 1/ \J x 2 + y 2 + z 2 + l 2 , and /3 is a constant. The form of Eq. Q21(l is similar to the equation 
found when Rosenbluth potentials [2(J are used in a gravitational potential <&. Note that the collision term has mixed 
velocity derivatives and that the resulting coefficient matrix from the discretization is non symmetric. We used central 
difference and a twenty- five point molecule to perform the discretization of Eq. Il21[| . see also Table HJ Here, we find 
the solution of the above equation in a six dimensional 'box' of length 1.22 units. At the borders we used a Dirichlct 
boundary condition of the form / = cxp(— x 2 — y 2 — z 2 ) exp(— v 2 — v 2 — v 2 ). 

We first started with coarse grid of four points per variable that leads to a matrix of 4096 2 elements (in this 
number we are not considering the border grid points given by the boundary condition) and (3=1. We found for this 
case that our method spent approximately 0.2 seconds to find the solution (all the calculations were perform with a 
Pentium IV of 1.8 GHz, the Fortran code was compiled with the Linux free compiler). We stop the iteration when 
^™ =1 (A* +1 ) < 10 -10 . The gmres routine and the public GMRES code [3(j spent approximately 0.35 and 3 seconds 
respectively to solve the system of equations with the same stop criteria, but it could be more if we choose wrong the 
restarting parameter. Here we are considering only the time spent to solve the coefficient matrix and not the time due 
to create the coefficient matrix and upload it into the code. In our code, we only upload 25 vectors of approximately 
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FIG. 4: Efficiency of our modified Stone method when compared to the gmres routine of Matiab that handies sparse matrices. 
When we increase the number of grid points our method becomes more efficient. To handle the 9 6 and 10 6 grid in our computer, 
we lower the value of the restarting parameter of the GMRES routine to the maximum possible value that avoids computer 
memory problems. 

For a grid of five points per variable some code compilers can not allow the storage of the coefficient matrix because 
is too large 5 12 , approximately 244 millions of elements. This was the case of the public code because it can not 
managed sparse matrices. For this number of grid points the modified Stone method was almost 7 times faster than 
the gmres routine. In Fig. ^we present the efficiency of our method compare to the gmres routine. Note that when 
we increase the number of grid points the efficiency also increases. For a grid of nine points per variable the gmres 
routine can not managed to find the solution in one iteration because we have computer storage problems. To handle 
this, we lower the restarting parameter to the maximum possible value that avoids this problem. For a grid of ten 
point per variable the number of elements in the coefficient matrix increases to 10 12 . Our code can managed this huge 
amount of data in a faster and efficient way. The time spent for this case was approximately 85 seconds. This should 
be the time spent for each time step in a non stationary problem, which is a great result considering the number of 
nodes in the grid and the precision attained with the stop criteria. The variation of a between the accepted limits 
< a < 1 slightly alter the time spent in finding the solution. In Fig. \5\ we present two graphs of the distribution 
function from the numerical solution of Eq. (|21|l for a grid of ten grid per variable. One is the surface distribution 
function on the plane x — y and the other is the velocity distribution function on the plane v y — v z ; they are plotted at 
different points of the grid. Note that these figures have not Maxwellian form because Eq. I|21[l is not in the collisional 
regime for the parameters considered in the collision term. This example was chosen only for didactic reasons. 

Note that when we used central difference for discretization of Eq. (|21|l . the only contribution to the main diagonal 
comes from the last right hand term, i.e. the term with j3. For some values of j3 < 1 our code diverge because 
our resulting coefficient matrix is not diagonal dominant. For the same value of /3, the gmres routine converges. 
The break in convergence at these values of (3 coincide with the appearance of negatives values on the distribution 
function solution found by gmres. We know that one of condition to attained a physical solution of the Fokker-Planck 
equation is that the distribution function has to be always positive. It is remarkable that for Eq. I|21|) our code 
diverge when physical solutions are not possible, this could be an indication that we are applying a wrong scheme for 
discretization or that we are not describing well the physical phenomena considered. As we said in the Introduction, 
other discretization scheme may be applied to avoid this problem. Thus, in our method the convergence is conditioned 
to the form of the functions present in and @, i.e. physical considerations; and to the difference scheme applied 
for discretization, i.e. numerical implementation. For a non stationary scheme the Crank-Nicolson decomposition 
in time is suggested because it has implicit character and it is usually more stable than other methods, but strictly 
speaking, the stability of the system has to be studied for each particular case considered. 



V. 



DOMAINS WITH CURVED BOUNDARIES 



The LU code was tested in the previous section with a structure-orthogonal grid but this does not mean that it 
can not be applied to more general geometries. To handle a non square domain we proceed as follows. First we make 
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FIG. 5: Numerical solutions of the distribution function of Eq. 12 1 1 . Top: surface distribution for the case with v x « 0.38, 
v y ~ —0.5, v z m —0.5 and z ~ 0.16. Bottom: velocity distribution for the case with x ~ —0.39, y ~ 0.06, z ~ —0.17 and 
v x « 0.28. 

a square structured grid with N no d e = ninjnknin m n n nodes, then we label them according to the index p denned in 
Later, the nodes that laid outside the boundary regions are not considered for the calculations. Now, lets see this 
procedure in an example. For didactic purpose we used only the one dimensional case of the Fokker-Planck equation 
in which we have two variables [x, v x ). In Fig. |H|is depicted in the two dimensional plane (x,v x ) a domain region f2 
with boundary Qg. The domain f2 is filled with a square structured grid that leads three classes of grid points: the 
interior points in which the normal discretization procedure can be done; the boundary points in which special care 
have to be taken when Dirichlet, Neumann or mixed boundary conditions are applied; and the exterior points that 
have to be neglected for the calculation, see Fig. |S| 

The incorporation of the Dirichlet boundary conditions using central differences for the boundary points of region 
f2 can be done using a Taylor expansion around the nearest nodes. For example, in Fig. f7| we take the nearest 
boundary points for point 1: the internal node (node 2) and the point at the boundary fig (point E); and make two 
Taylor expansions around point 1. These expansions give us (our conventions are: partial derivative with respect to 
the coordinate x denoted by (, x); Ax is the discretization interval in the x direction), 



12 




FIG. 6: Application of the incomplete LU decomposition to a general domain Q with curved boundaries fin- The square 
structured grid leads three classes of grid points: interior points (full circles), boundary points (empty circles) and exterior 
points (triangles). 





FIG. 7: Schematic representation of a five-point molecule for the finite difference approximation of the derivatives in the 
plane (x,v x ) for a boundary point. Below, we see different boundary situations that can occur in the numerical computation 
depending on the discretization grid and boundary S. 



(bAx) 2 

= *i + bAxV, x + . 

(Ax) 2 

* 2 = - Ax^. x + ^—J-^ xx 



: + <3[Ax 3 ], 
0[Ax 3 }, 



(22) 



eliminating the second order derivative between these equations we can express the partial derivative at a boundary 
point as, 



Ax 



1 



6(1 + 6) 



b 1-6 

*2 



1 + 6 



C>[Aa 



(23) 



When b = 1 we recover the usual expression for the central derivative. The same procedure can be applied to the 
northern boundary, between node 3 and point A, and to the different types of boundary nodes in Fig. [7| Also, partial 
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derivative of higher orders and mixed derivatives can be found in a similar way. Furthermore, it is also possible to 
implement Neumann boundary conditions in an irregular boundary using finite difference, see for instance |3l|. 

The exterior points of Fig. are strictly necessary to maintain the ordering of the nodes in a domain with curved 
boundaries. This ordering is needed by the code to operate normally but they do not enter into the calculations. To 
implement this condition, we have to set the values of the elements £pn> ^P' ^pf]> ^ P an< ^ ^ P °^ ^ ne external 
grid point (each exterior point has a position p from the one dimensional storage index) equal to zero. 

An application of the LU method in a two dimensional Fokker-Planck equation (x, y, v x , v y ) with curved boundaries, 
as well as the incomplete LU decomposition for this case can be seen in [2{|. Here, the distribution function of a 
stationary gravitational thin disk is calculated. Note that in order to obtain the elements of the LU decomposition 
for the two dimensional case, we have to start the calculations from the beginning, i.e. we can not used the elements 
found in the three dimensional case. 



VI. CONCLUDING REMARKS 



We have developed a variation of the incomplete LU decomposition proposed by Stone that solves the three 
dimensional linear Fokker-Planck equation. The method presented can manage the large set of linear equations that 
appears from the discretization procedure. The convergence of the iterative process is done in a fast and effective 
way. Also, this method can be easily adapted to support irregular boundaries with Dirichlet, Neumann or mixed 
boundary conditions and can be used to follow the evolution of a distribution function for non stationary equations. 
In this case a Crank-Nicolson discretization in time is recommended because it has implicit character and is more 
stable than other methods, but strictly speaking, the stability of the system has to be studied in each particular case 
considered. The good properties that our method shares and the lack of methods that handle the large amount of 
data (given by the Fokker-Planck equation with six independent variables) make the method presented in this article 
worthy and advantageous. In general, the algorithm presented in the article can solve the system of equations that 
arises from the discretization of the equation or similar equations (with or without the mixed derivatives), but we 
have to keep in mind that its convergence is conditioned to the form of the functions v, A(x., v), £?(x, v), C(x, v), 
D(x, v) and Eij(x., v) that appear in the collision term @. 



APPENDIX A: ELEMENTS OF M = LU AT NODE p 



In this section we present how the elements of the matrix M are related to the elements of the upper U and lower 
L matrices. 



1V1 DD3 — ^DD3i 

l\fP _ TP j]-p-{nijklm+nijkl) 

1V1 DD3\N ~ DD3 N > 

MP — TP T T P-( n i3ki m +n ijk i) 

1V1 DD3\E ~ DD3 E > 

MP — TP T TP-( n i]ki m +n,jki) 

1U DD3\T ~ DD3 T ' 

MP — TP jjP-( n ijki m +rtijki) 

1V1 DD3\UI ~ lj DD3 u Ul > 

MP —TP -L TP TTP-{nijkim.+nijki) 
1V1 DB3 ~ DB3 lj DD3 u UB2 ' 

ijp TP TTP~( n ijklm+nijk) 

M DB3\N ~ DBS N ' 

1,/rp _ T p jjP-(nijklm+nijk) 

DB3\E ~ DB3 E > 

M p _ T p TT P~(n ijk i m +n ijk ) 

1U DB3\T ~ DB3 T ' 

nr p _ jP i TP jjP-(nijkl m +nijk) , T p T jP~(nijkl m +nijki) 
D3 D3 ' DB3 Ul ' ^DD3 U U2 > 

M p —TP jjP-nijkim 
IVI D3\N ~ n D3 u N ' 

i\/rP _ TP TjP- n ijki m 
IVI D3\E ~ lj D3 u E 1 

MP —TP TTP~ n i3kl m 

1U D3\T ~ n D3 u T > 

nrP —TP J- TP TjP- n ijki m . T p TjP-{nijki m +nijki) 
DT3 DT3 D3 Ul ' lj DD3 u UT2 ) 
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M p _ T v TTP-(nijklm-riijk) 

DT3\N ~ DT3 N ' 

tljP _ jp TjP-{nijklm-nijk) 

IV1 DT3\E — n DT3 u E i 

M p _ jp T jP-{riijkl m -n ijk ) 

1V1 DT3\T ~ 1J DT3 U T ' 

M p _ jp jjP-(ni jk i m -n ijk ) 

1V1 DT3\U1 ~ n DT3 u Ul ' 

M p _ jp T TP-(n ijklm +n ijk ) 

1V1 DB3\UB2 ~ ^DB3 U UB2 ' 

M p _ T p TJ p-n ijklm T p TT p-(n ijk i m +n ijk ) 

1V1 D3\UB2 ~ lj D3 u UB2 "r lj DB3 u U2 ' 

M p — TP I rP T jP-{nij k i m -n ijk ) jV TJ p-n ijk i m jV j jp-{n ijk i m +n ijk ) 
1V1 DU3 ~ ^DU3 ' ^DT3 U UB2 ' t1j D3 u U2 * 1j DB3 u UT2 ' 

M P _ jP TjP-(n ijklm -n ijkl ) 

1V1 DU3\N ~ n DU3 u N ' 



M p _ jp TjP-(jiij klm -n ijkl ) 

1V1 DU3\E ~ ±J DU3 U E ' 

i/P _ TP TjP-(jiij klm -n ijkl ) 

1V1 DU3\T ~ 1J DU3 U T ' 

M p _ jp TTP-(n ijklm -n ijk i) jp T ]-p-(n ijk i m -n ijk ) jV T j-p-n ijk i m 

1V1 DU3\U1~ DU3 U U1 ' ^DT3 U U2 ' ^ D3 U UT2 ' 

l/P _ TP TjP-(ji ijklrn -n ijk ) 

1V1 DT3\UT2 ~ ^DT3 U UT2 ' 

Tiyrp _ T p T TP~(nij k l m ~n ijk i) 

1V1 DU3\UB2 ~ n DU3 u UB2 > 

Tiyrp _ T p TTP-(ni jk i m -n ijkl ) 

1V1 DU3\U2 ~ lj DU3 u U2 ' 

Tiyrp _ T p TTP-(ni jk i m -n ijkl ) 

1V1 DU3\UT2 ~ ±J DU3 U UT2 > 

tit-p _ T p jjP-(n ijk im+nijkl) 

1V1 DD3\UD3 ~ ^DD3 U UD3 > 

titP J V I TP rrP-("ijfclm+"ijfc) . jp TjP-inijklm+Tlijkl) 

DB2 DB2 ' DB3 UD3 ' ^DD3 U UB3 ' 

nfP _ TP TjP-{ n ijkl+n iik ) 

DB2\N ~ 1J DB2 U N ' 

n/rp _ jp T T p-(nijkl+n ijk ) 

1V1 DB2\E ~ ^DB2 U E ' 

l/P _ tP jjP-i^ijki+nijk) 

1V1 DB2\T ~ ^DB2 U T ' 

M p _ jP , jP T jP-{nij k i+n ijk ) T p TJ p-n ijklm T p TJ p-(n ijklm +n ijkl ) 

1V1 D2 — ^D2 ' lj DB2 u Ul ' ^ D3 U U D3 ' ^DD3 u U3 ' 

mp — jp TiP- n im 

1V1 D2\N ~ n D2 u N ' 

M P _ TP TjP-Hijkl 

1V1 D2\E ~ lj D2 u E ' 

M p —TP TjP~ n i3kl 

1V1 D2\T ~ n D2 u T ' 

titP _ T p i TP jjP- n i3kl | tP TjP-(nij k i m -n ijk ) T p T Tp-(n ijk i m +n ijk i) 
1V1 DT2 DT2 ' D2 U U1 ' ^DT3 U UD3 ' lj DD3 u UT3 ' 

M p —TP TjP-(nijkl-nijk) 

1V1 DT2\N ~ ^DT2 U N ' 

l/P _ tP jjP-{nij k i-n ijk ) 

1V1 DT2\E ~ lj DT2 u E ' 

M p _ jp T jP-{n ijk i-nij k ) 

1V1 DT2\T ~ DT2 T ' 

l/P _ tP TTP-(nij k i~n ijk ) 

1V1 DT2\UI ~ n DT2 u Ul ' 

i/P _ tP TjP-(ni jk i+n ljk ) T p jjp-inijkim+riijk) 

1V1 DB2\UB2 ~ n DB2 u UB2 ' ^DB3 U UB3 ' 



i,fP — TP i TP TjP- n ijkl i jp TTP-(ni jk i+n ijk ) T p jjp-nij klm jV r Tp-(n ijk i m +n ijk ) 

1V1 D1 ~ n Dl ' lj D2 u UB2 ' ^DB2 U U2 ' ^D3 U UB3 ' ^DB3 U U3 

M p —TP T jP- n ijk 
IVI D1\N ~ n Dl u N ' 

MP — TP TTP~ n Hk 
M D1\E - L D\ V E 

MP — T p TjP~ n Hk 
1V1 DI\T ~ ^Dl u T ' 
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M P ]N = L p B U p N ~ n ^ 

M p BlE = L* B Ul- ni \ 
M P W = L P W , 

1V1 W\N ~ W N > 

M p = L% 

]\/rP — TP \ TPttP- 1 i TP JjP^ n i i TP TjP- n ii i TP jjP~ n iik . T p TjP-{n ijk i-n ijk ) jV TT p-n ijkl 
P P ' S N +^W U E "+" Ll B v T +lj Dl u Ul + ^DT2 U UB2 + ^ D2 U U2 

,tP T jP-(nijki+nij k ) , T p T TP-(nij k i m -n ijk i) T p j T p-(n ijk i m -n ijk ) ,jp jjP-n ijk i m 
' lj DB2 u UT2 "T ^DU3 U UD3 ' ^DT3 U UB3 ~ rlj D3 u U3 

I TP TTP—{ n ijklm+Tlijk) I TP T jP~ ( n ijklm+Tlij k l) 

<^DB3 U UT3 lj DD3 u UU3 ' 

M P N = L p p U p Nl 

M p = L P P U P , 

M p —T p n p ~ nj 

W\T ~ W T ' 

m s\t = l s^t \ 
MP = L p p U p , 

MP -T p TT v ~ n%3 
M B\U1 - L B U U1 

MP -TP TT p ~ n ' 
1V1 W\U1 ~ n W u Ul > 

M S\U1 = ^s^c/i 1 ' 

MP — TP TTP i TP TjP-{nijkl-riij k ) yp jjp-n iik i jp T jp-(nij k i m -n ijk ) jp jjp-n iiklm 
1V1 U1 ~ 1J P u Ul " r ^DT2 U U2 "r" ^D2 U UT2 ~^ ±J DT3 U U3 ' Ju D3 u UT3 ' 

_ T p T jP-(nij k i-n ljk ) jp jjp-{riij klm -n ijk ) 
1V1 DT2\UT2 ~ ^ DT2 U UT2 ' lj DT3 u UT3 ' 

M P — T p TT p ~ n ^ k 

1V1 DI\UB2 ~ lj D\ u VB2 ' 

M P — T p TT p ~ n%j 

1V1 B\UB2 ~ lj B u UB2 ' 

M P — T p TT p ~ ni 

1V1 W\UB2 ~ lj W u UB2 > 

M p - T P TI P ~ 1 
1V1 S\UB2 ~ n S u UB2> 

Tiyrp _ tPtjP i TP jjP-riijk , jp TTP-(nij k i m -n ijk i) jV TT p-(n ijklm +n ijk ) 
1V1 UB2 — ^P U UB2 ' L 'UV J U2 ' rJ ^DU3 u UB3 ' lj DB3 u UU3 ' 

MP -T v TT p ~ n * 3 
1V1 B\U2 ~ ^B U U2 ' 

M P -TP jjP- n i 

1V1 W\U2 ~ n W u U2 ' 

M p - tptt'p- 1 

1V1 S\U2 ~ lj S u U2 ' 

MP — TP TTP i TP TjP- n ijk , jp T TP-{ni jk i m -n ijkl ) T p TT p-n ijklm 
1V1 U2 ~ lj P u U2 ' ±J D1 U UT2 ' ^DU3 U U3 ' rlj D3 u UU3 ' 

MP -T v TT p ~ n%j 

1V1 B\UT2 ~ ^B U UT2 ' 

M P -TP jrP- n i 

1V1 W\UT2 ~ ±J W U UT2 ' 

M p - T V TTP~ X 

1V1 S\UT2 ~ lj S u UT2' 

iiyrp _ tPtjP i tP TTP-(ni jk i m -n ijkl ) jP TT p-(n ijklm -n ijk ) 
1V1 UT2 ~ lj P u UT2 n DU3 u UT3 ' ^ DT3 U UU3 ' 

MP — TP TjP-{nijki m -nij k i) 

DU3\UU3 ~ 1J DU3 U UU3 ' 

MP — TP jjP-(jiijki+nij k ) 

1V1 DB2\UD3 ~ lj DB2 u UD3 ' 

MP —T p TTP~ n ii>*l 

1V1 D2\UD3 ~ n D2 u UD3 ' 

MP — TP TjP-{nij k i-n ijk ) 

1V1 DT2\UD3 ~ n DT2 u UD3 ' 

MP — TP TT P-(nijkl+n.ijk) 

1V1 DB2\UB3 ~ ^DB2 U UB3 ' 



10 



MP — T p TT p ~ n ^ k _L TP TT p ~ ni i kl J_ TP j]-p-( n ^ki+n ijk ) 

IV1 D1\UD3 ~ L 'D\ tJ VD3 ' n D2 u UB3 ' lj DB2 u U3 > 

M B\UD3 ~ ^B U UD3 ' 

MP -T p Tl v ~ rii 

1V1 W\UD3 ~ lj W u UD3 ' 

M p — T V TT V ~ X 

S\UD3 ~ ^S U UD3' 

MP — TP TIP I fP rrP-("«W-«y*) , rP jjP~ri ijk i T p TJ p-{n ijkl +n ijk ) 
1V1 UD3 ~ lj P u UD3 t n DT2 u UB3 ' ^D2 U U3 ' lj DB2 u UT3 > 

»,fP _ rP Tj-p-(ni jk i-n ijk ) jV TJ -p-n ijk i 

1V1 DT2\U3 ~ n DT2 u U3 ' ^ D2 U UT3 ' 

MP — T p Tj-P-( n ijki-riij k ) 

1V1 DT2\UT3 ~ ^DT2 U UT3 ' 

M p - T p TT p ~ ni:jk 

1V1 DI\UB3 ~ lj Dl u UB3 ' 

MP - T p TJ v ~~ n% i 

1V1 B\UB3 ~ lj B u UB3 ' 

MP -T p TT p ~ nj 

W\UB3 ~ ^W U UB3 ' 

M p — TPTJP^ 1 

1V1 S\UB3 ~ n S u UB3' 



p—nij k T p jjP—(nij k i+ni jk ) 



MP — T v TIP 4- T v TT V ~ ijk 4- TP Tl p ~ [ 
1V1 UB3 ~ n P u UB3~ t n Dl u U3 ~ rIj DB2 u UU3 

M p - T p TT p ~ ni] 
1V1 B\U3 ~~ ^B U U3 ' 

MP —T p TT p ~ nj 
W\U3 ~ ^W U U3 ' 

M P - T P TT P ~ 1 
1V1 S\U3 ~ lj S u U3 ' 

MP — T P TT P J- TP TT p ~ ni i k J_ TP TT p ~ nt i k > 

IV1 U3 ~ n P u U3 "r" n Dl u VT3 ' lj D2 u UU3 ' 

M p - T p T7 p ~ nij 

1V1 B\UT3 ~ n B u UT3 ' 

M P — TP TjP- n i 

1V1 W\UT3 ~ n W u UT3 ' 

M P — rPrrP- 1 
1V1 S\UT3 ~ 1j S u UT3' 

MP — TPTTP _l TP TTP-( n ijki-n ijk ) 
1V1 UT3 ~ lj P u UT3 n DT2 u UU3 ' 

M P — T p TJ p ~ ni3k 

1V1 DI\UU3 ~ 1J Dl tJ UU3 ' 

MP - T p TTP^ 11 ^ 

1V1 B\UU3 ~ lj B u UU3 ' 

MP - T p T7 p ~ nj 

1V1 W\UU3 ~ ^W U UU3 ' 

yrP - TPTTP" 1 

1V1 S\UU3 ~ ^S U UU3> 



Kvz = VpUlvs- ( A1 ) 



APPENDIX B: EXPLICIT FORM OF EQ. © 



In this section we present the explicit form of w for each grid node. 



OdD3^ DD3 + OdB3^ DB3 + Od3^D3 + O DT3^ DT3 + OdU3^ DU3 + OdB2^ DB2 + Od2^ D2 + OdT2^ DT2 
+ D1 ^ D1 + B ^B + Ow^W + Os^S + Opfp + N ^ N + Oe^E + Ot*T + + UB 2^UB2 

+Ou2^U2 + OuT2^UT2 + OjjDZ^UDZ + OuB3^UB3 + Ou3^U3 + OuT3^UT3 + OijU3^UU3 + MdD3\N^ DD3\N 
+M DD3 \ E 'i' DD3\E + ^DD3\T^ DD3\T + -^DD3|yi*DD3|C/l + M D B3\ N ^ DB3\ N + ^DB3\E^ DB3\E 
+M DB3 \ T ^ EB3\T + MD3\N^ D3\N + Md3\E^ D3\E + Md3\T^ D3\T + ^DT3\ DT3\N + ^DT3\E^ DT3\E 
+M DT3 \ T iS> DT3\T + ^DT3\U1^ ' DT3\U1 + ^DB3\U B2^ DB3\U B2 + ^DB3\U2^ DB3\U2 + ^ D3\U B2^ D3\U B2 
+M DU3 \ [N i l! EU3\N + M DU3 \ E ^ DU3\E + ^DU3\T^ DU3\T + ^D3\UT2^ D3\UT2 + M DT3 \ U2 ^ DT3\U2 
+M DU3 \m^ DU3\U1 + M DT3 \ UT 2^ DT3\UT2 + M D U3\UB2^ DU3\UB2 + M D U3 1 U2 ^ DU3\ U2 + -^D(73|(7T2^'D(73|(7T2 
+M DD3 \ UD3 ' i i> DD3\UD3 + M DB2 \N^DB2\N + M DB2 \E^ DB2\E + M DB2 \ T ^ ' UB2\T + M D2 \ N ^ D2\ N + M D2 \ E ^ D2\E 
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+M E2 \T^ D2\T + ^DT2\N^ DT2\N + MdT2\E^ DT2\E + MjjT2\T^ DT2\T + MbT2\UI^ DT2\U\ 

J rM DB 3\UB3^ DB3\UB3 + M DB 2\UB2^ DB2\UB2 + ^Dl\N^ Dl\N + Mdi\E^D1\E + Moi|T* ' Dl\T + M B \N^B\N 
+Mb\E^ B\E + M\V\N^>W\N + M S \E^S\E + M W \ T ^ W \ T + M S | T \1> S | T + M B \U1^B\U1 + 
+Ms\ui^S\Ul + M DT2 \UT2^ DT2\UT2 + M DT3 \ UT3 ^ DT3\UT3 + M- D1\U B2^ ' D1\U B2 + M B \UB2^ B\VB2 
+M W \UB2^W\UB2 + M s\u B2^ S\U B2 + M B \U2^ B\U2 + M W \ U2 1 $> W \ U2 + ^S\U2^S\U2 + M B \ UT2^ B\UT2 
+M W \UT2^W\UT2 + M s \UT2^ S\UT2 + M DU3 \ UU3 ^> ' DU3\UU3 + M DB 2\UD3^ DB2\UD3 + M d2 \jj d3 *5> D2\UD3 
+M DT2 \ um '$> DT2\UD3 + M DB 2\UB3^ DB2\UB3 + ^ D1\U D3^ D1\U D3 + M D2 \UB3 ^D2\UB3 + M DB 2\U3^ DB2\U3 
+Mb\UD3^ B\UD3 + ^W\UD3^W\UD3 + M S\U D3^ S\U D3 + M DT2 \ U3 *S> £>T2|(73 + -^D2|(7T3 1 I'D2|i7T3 
+M DT2 \ UT3 ^ DT2\UT3 + -^D1|(7B3 V I / D1|(7B3 + -^B|(7,B3 , I'b|(7B3 + ^W\U B3^W\U B3 + Ms|(7B3 V I , S|E/.B3 
+M B \U3^ B\U3 + | (73 ^1^1(73 + Ms\U3^ S\U3 + J ^B|(7T3 V I / B|(7T3 + -^W \UT3^W\UT3 + -^S|(7T3^S|(7T3 

+-^D1|C/(73 V I'D1|(7(73 + M B \UU3^B\UU3 + -^W|(7(73^ , W|(7(73 + -^51(7(73^51(7(73 ~ 0. (Bl) 

APPENDIX C: O AS A LINEAR COMBINATION OF M 

In this section we present the elements of the matrix O as a linear combination of the elements of matrix M using 
the relations (|10|) . 

Odz>3 = — u(M DD3 \ N + M DD3 \ E + M DD3 \ T + M DD3 \ui + M DD3 \ UD3 ), 

OdB3 — ~0i(M DB 3\N + M DB 3\E + MdB3\T + MdB3\UB2 + M DB 3\U2 + M DB 3\UB3): 
Od3 — —a(Mo3\N + MD3\E + MD3\T + M D3 \ UB 2 + M D3 \ UT2 ), 

Odt3 = —oi{M DT3 \ N + M DT3 \ E + M DT3 \ T + M DT3 \ui + M DT3 \ U2 + M DT3 \ UT2 + M DT3 \ UT3 ), 

OdU3 = — a {M DU3 \ N + M DU3 \ E + M DU3 \ T + Mbu^u! + M DU3 \ UB2 + M DU3 \ U2 + M DU3 \ UT2 + M DU3 \ UU3 ), 
OdB2 = — a (^DB2|7V + MdB2\E + MdB2\T + MdB2\UB2 + MdB2\UD3 + MdB2\UB3 + MdB2\U3)i 
Od2 = —Ot(Mij2\N + Md2\E + Md2\T + ^D2\UD3 + Md2\UB3 + M D2 \ UT3 ), 

OdT2 = ~ a {MDT2\N + -^_DT2|B + MdT2\T + ^7JT2|(71 + ^DT2\UT2 + M DT2 \ UD3 + M DT2 \ U3 + M DT2 \ UT3 ), 
Odi = —a(M Dl \ N + M Dl \ E + M D i\ T + Md\\UB2 + Mdi\ud3 + M m \UB3 + Mdi\uu3), 

Ob = —a(M B \N + Mb\E + Mb\U1 + M B \UB2 + M B \JJ2 + Mb\UT2 + Mb\UD3 + Mb\UB3 + Mb\U3 + -^B|!7T3 
+M B \UU3)i 

Ow = — ot{M w \ N + M W \ T + M W \ui + M w \ub2 + M W \ U2 + M W \ UT2 + M W \ UD3 + M W \ UB 3 + M W \ U3 + M W \ UT3 
+M W \ UU3 ), 

Os = —a(Ms\E + Ms\T + Ms\Ul + Ms\UB2 + Ms\U2 + Ms\UT2 + M S \ UD3 + M S \UB3 + Ms\U3 + M S \UT3 
+M S \UU3), 

Op = Ot{M DD3 \ N + MdD3\E + Mdd3\t + Mdd3\ui + Mdb3\n + Mdb3\e + Mdb3\t + Md3\n + Md3\e + M D3 \ T 
+M DT3 \ N + M DT3 \ E + M DT3 \ T + M DT3 \ [U1 + M DB3 \ub2 + M DB 3\U2 + M D3 \ UB 2 + M DU3 \ N + M DU3 \ E 
+Mdu3\t + Md3\ut2 + M DT3 \ U2 + M DU3 \ui + M DT3 \ UT2 + M DU3 \ UB2 + M DU3 \ U2 + M DU3 \ UT2 

+M E D3\UD3 + Mbb2\N + MdB2\E + M E B2\T + Md2\N + Md2\E + ^D2\T + M E T2\N + MdT2\E + MbT2\T 
+M DT2 \ui + M DB3 \UB3 + MDB2\UB2 + M m \N + M m \E + M D \\T + M B\N + Mb\E + M\V\N + M S \ E 

+M W \ T + M S \ T + M B \ui + M\v\ui + Ms\ui + M DT2 \ UT2 + M DT3 \ UT3 + M D1 \ UB2 + M B \ub2 + ^w\UB2 

+M S \ UB2 + M B \U2 + M\Y\U2 + M.S\U2 + M B \UT2 + M W \ UT2 + M S \ UT2 + M DU3 \ UU3 + M DB2 \UD3 
+M D2 \ UD3 + M DT2 \ UD3 + M DB2 \UB3 + -^B1|(7B3 + MD2\UB3 + M DB 2\U3 + ^B\UD3 + M\V\UD3 + M S \ UD3 

+M DT2 \ U3 + M D2 \ UT3 + M DT2 \ UT3 + M D i\ UB3 + M B \ub3 + M\v\UB3 + Ms\UB3 + M B \U3 + M W \ U3 

+M,s\U3 + M B \UT3 + M W \ UT3 + M S \ UT3 + M D i\ UU3 + M B \UU3 + M\V\UU3 + ^S\UU3)i 

On = -a(M DD3 \ N + M DB3 \n + m D3\n + m dt3\n + Mdu3\n + ^db2\n + m D2\n + m dt2\n + Mdi\n + M b \n 
+M W \ N ), 

O e = —a(M DD3 \ E + M DB3 \ E + M D3 \ E + M DT3 \ E + M DU3 \ E + M DB 2\e + M D2 \ E + M DT2 \ E + M D1 \ E + M B \ E 
+M S \ B ), 
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Or — —ot{M Dm \ T + M DB3 \ T + M m \ T + M DT3 \ T + M DU3 \ T + M DB2 \t + M D2 \t + M DT2 \t + M m \ T + M W \ T 
+M S \ T ), 

Oui = -a(Mr>D3\ui + M DT3 \ V1 + M^^m + M DT2 \ui + + M W \ U1 + Ms\ui), 

OuB2 = —Oi(M DB3 \ UB2 + M D3 \ UB2 + M DU3 \ UB2 + M DB2 \ UB2 + M D i\ UB2 + M B \UB2 + M\V\UB2 + Ms\UB2), 

Ou2 — —oc{M DB3 \ U2 + M DT3 \ U2 + M DU3 \ U2 + M B \ U2 + M W \ U2 + M S \u2)i 

Out2 = —ol{M D3 \ UT2 + M DT3 \ UT2 + M DU3 \ UT2 + M DT2 \ UT2 + M B \ VT2 + M W \ UT2 + M S \ UT2 ), 

OuD3 — —a:{M DD3 \ UD3 + M DB2 \ UD3 + M D2 \ UD3 + M DT2 \ UD3 + M Dl \ UD3 + M B \ UD3 + M W \ UD3 + M S \ UD3 ), 

Oub3 = —a(M DB3 \ UB3 + M DB2 \ UB3 + M D2 \ UB3 + M D i\ UB3 + M B \ UB3 + M W \ UB3 + M S \ UB3 ), 

Ou3 = -a{M DB2 \ U3 + M DT2 \ U3 + M B \ U3 + M W \ U3 + M S \u3), 

Out3 = —oc{M DT3 \ UT3 + M D2 \ UT3 + M DT2 \ UT3 + M B \ UT3 + M w \ut3 + M S \ UT3 ), 

Ouu3 — -ol(M DU3 \ UU3 + M D i\ UU3 + M B \ UU3 + M W \ UU3 + M S \ UU3 ). (CI) 



APPENDIX D: EXPLICIT FORM OF THE FUNCTIONS C [x] AND K [x] 
In this section we show the functions Cpf] and K[x] that appear in the final expression (|11|1 of the LU decomposition. 

_ TP TTP—{ n ijklm+nijkl) 

^DB3 — L DD3 V VB2 > 

si jp jjP-inijkim+nijki) T p TT p-(n ijklm +n ijk ) 

L-D3 — ^DD3 U U2 +lj DB3 U Ul ' 

_ TP rrP _ ( n *3*l"»+ n *3*l) I TP TTP~ n ii>'lm 

ODT3 — ^DD3 U UT2 + L D3 U U1 > 

n t p Tj-p-(n ijk i m +n ijk ) T p TJ -p-n ijk i m T p j T p-(n ijklm -n ijk ) 

^DU3 — ^DB3 U UT2 ' rlj D3 u U2 ' ^DT3 U UB2 ' 

(-1 _ T p TT P-(nij k i m +n ijk i) T p TT p-(n ijk i m +n ijk ) 

^DB2 — ^DD3 U UB3 + ^DB3 U UD3 ' 

r , _ r p TT P-(nij k i m +n ijkl ) T p TT p-n ijklm T p TT P~{n ijk ,+n ijk ) 

^D2 — ^ DD3 U U3 -I- ^D3 U UD3 ^^DB2 U U1 ' 

n _ T p TT p-(ni jk i m +n ijkl ) T p TT p-(ni jklm -n ijk ) jV TT p-n ijkl 

^DT2 — ^DD3 U UT3 ^DT3 U UD3 ^D2 U U1 ' 

(~i _ jp T]-p-{nij k i m +n ijk ) T p TJ -p~n ijk i m jP j-j-p-(n ijkl +n ijk ) jV TJ -p~n ijkl 

^Dl — lj DB3 u U3 ' ^D3 U UB3 ~^^DB2 U U2 ' ^ D2 U U B2 i 

f-i _ T p TTP-(nijkl-n ijk ) T p j-jP-n ijkl jV r j-p-(n ijklrIl --n ijk ) jV jjp-n ijklm 

^Ul — lj DT2 u U2 "T n D2 u UT2 ~ tl ^DT3 u U3 ' ^ D3 U UT3 ' 

n _ T p TT P~n ijk T p T] -p-(n ijk i m ~n ijk i) jV T r p- (n ijklm +n ijk ) 

^UB2—^ D iU U2 +^DU3 U UB3 +^DB3 U UU3 ' 

— TP T[P— n ijk , jP Tj-p—(nijklm—nij k i) jp jjP—riij klm 
^V2—^UV J UT2 ' ± - u DU3 u U3 <^D3 U UU3 ' 

(~i _ T p TT P-{nij klm -n ijkl ) T p TT P~(n ijklm -n ijk ) 
lsUT2 — ^DU3 U UT3 ^ DT3 U UU3 ' 

n t p T jP-(n ijk i-n ijk ) T p jjp-nij k i jV TT p-(n ijkl +n ijk ) 

<^UD3 — ^DT2 U UB3 ' lj D2 u U3 ^ lj DB2 u UT3 ' 

n _ r p Tr p-nijk | T p TT p-(nij k i+n ijk ) 

^UB3—^D1 U U3 + lj DB2 U UU3 ' 

si _ TP jrP-nijk , T p jjP-riijki 

^U3 — 1J Dl u UT3 ^ 1J DT J UU3 ' 

n _ T p TT p-(nij k i-n ijk ) 

^UT3 - h DT2 U UU3 ' 

Cdd3 = C B = Cw — Cs = Cn = Ce = Ct = Cuu3 — 0, 

TT p-(nij k i m +n ijk i) TT p-(n ijk i m +n ijk i) TT p~(n ijk i m +n ijk i) TT p-(n ijk i m +n ijk i) TT p~(n ijk i m +nii k i) 
J^DD3 = U N + U E +U T + U U1 + U UD3 ) 

is T TP—{nijklm+nijk) | TTP—(nijklm+riijk) . ttP— (nijfcim+nijfe) . rjP— ("ijfclm +nijk ) . tjP~- (nijfelm+iijfe) 

&DB3 = U N + U E + U T + U UB2 + U U2 

,T T P-(nij k i m +n ijk ) 

' U UB3 ' 
7^ _ TT p-n ijklm T Tp-n ijk i m jjp-riij klm Tr p-n ijk i m TT p—n ijk i m 
&D3-U N +U E +U T +U UB2 +UuT2 

TT p-(n ijk i m -n ijk ) TJ -p-(n ijk i m -n ijk ) Tr p-(n ijk i m -n ijk ) TT p-(n ijklm -n i:jk ) Tr p~(n ijk i m -n ijk ) 
J^DT3 = U N + U E +U T + U Ul + U U2 

,TTP-{ n ijklm-n ijk ) jjP-{nij k i m -ri iik ) 
< U UT2 "r U UT3 ' 

jjP—{nij k i m —ni jk i) TTp—(nij k i m —n ijk i) t Tp—(n i]k i m —riij k i) j-j-p—(nij k i m —nij k i) T j-p—(nij k i m ~nij k i) 

' N "r U E ' T U U1 + U UB2 

.TjP—{nijkim—riijkl) , ttP— {ntjhlm — "ijfci ) , j jP— ("ijfcim — "ijfel ) 
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tv- _ TT p-(nij k i+Hijk) | TTP-(nijki+nijk) . jjP— ("ijil +"ijfc ) . r rt>— (nyn+ny*) . TT p-(nijki+nijk) , j-t-P— (nyfci+ny*) 

-K.DB2 - (Vjv + + + ^1752 + ^£7.03 + ^[733 

, 7-rP— (nijki+riijk) 
+ U (73 ' 

tv- TT p—riijkl | T TP—riijkl | TjP—nijkl . TT p—nijkl . TT P~nijki jjP—riijkl 

A B2 = fjy + (7 E + <7 T + + + (7rj T3 , 

iv- _ TjP—fjiijki—nijk) | rrP— (niju— ny*) , r rP— ("ijfcl — nijfc ) i TTP—(riijki—nijk) . TT p—(nijki—nijk) , tjP— ("iifc! — ny*) 

«-DT2 - f ^ + (7 £ + (7 T + f ;71 + + (7^3 

,TTP—(nijki—nij k ) , jjp—ijiijki—riijk) 
+ u (73 ' U [7T3 ' 

„ TT p-n ijk TT p-n ijk TT p-n ijk TT p-n ijk TT p-n ijk TT p-n ijk TT p-n ijk 

&D1 -u N +u E +u T + u UB2 + u um + u um + , 

TV- _ TT P-n t j TT P-n tj jrP-nij TT p-nij , TT P-nij TT p-n tj TT p-n t] TT p-n t] jjP-nij rrP-ntj TT p-n t] 

K B - U N +U E + V ux + 67^2 + (V^ + (7^2 + (7^3 + (7^3 + (7^ + (7rj T3 + (7^3 , 

„ _ jjp-Tlj rrP-«j , T-rP-Wj , TjP-nj , rrP-Kj , TjP-nj , rrP-"j , rrP-Kj , fiP-lj , TrP-rij , rrP-"j 

J\\\r — u N ~r C T t~ Crjj t L^(7B2 "r u U2 ' U UT2 ' U UD3 ' U UB3 ' u U3 ' U UT3 + u UU3 ' 

K _ TJP-I , 7JP-1 , JJP-1 , TJP-1 , TJP-1 , TJP-1 , JJP-1 , 7JP-1 , 7JP-1 , JJP-1 , 7JP-1 

n S — v E T u T ^ (71 ^ 17 [7B2 ^ 17 C/2 ^ l7 !7T2 ' ^[7153 ^ U !7B3 ^ u U3 ^ u UT3 ' u UU3> 

_ T p TT p-{riij k im+nij k i) T p j-j-p-inijkim+riijk) , T p jjp-riijklm , r p TT p-(nijklm-nijk) 
J\N — ^ DD3 U N ' ^ DB3 N ' r - L D3 U N ' DT3 N 



1 rP T jP-(ji ijklm -n iikl ) T p TT p-(n ijk i+n ijk ) jP TT p-n ljkl jP TT p-(n ijkl -n ijk ) f p TJ p-n ijk 
~ t2j DU3 u N ~^^DB2 U N ~ rn D2 u N ' DT2 N ~ t -^Dl u N 



1 rP TjP-nii , rP rrP-"j 

TV- _ rP TjP—('n i j k i m +n i j k i) T p TT p—{n z j k i m +nij k ) f p TT p~nij k i m T p ttP— ("tjfcim— 

KE~^DD3 U E +1j DB3 U E + L D3 U B DT3 E 

_l_rP TjP-(nijkl m -n ijk i) . T p TT p-(n ijk i+ni jk ) jV TJ p-n ijk i jP T T p-{n ijk i -n ijk ) r p TJ -p-n ijk 
+lj DU3 U E +^DB2 U E "+" ^ D2 U E "+" ^ DT2 U E + ^ Dl u E 

7v- _ r p TT p-(nij k i m +n ijk i) T p TT p-(n z j k i m +n ijk ) T p TT p-n ijk i m T p TT p-(n ijk i m ~n ijk ) 
J±T — ^DD3 U T +lj DB3 U T +^D3 U T ' DT3 T 

1 rP jjP—{nijMm—nijkl) , j-p TTP-(ni jk i+n ijk ) jp jjp-n ijk i jp TrP-(iy«-B«l) . rP j-j-p-n ijk 
~ r±J DU3 u T "r lj DB2 u T ' ^ D2 U T ' ^ DT2 U T ' ^ Dl u T 

TV" rP TjP—(nijklm+nijkl) , 7-p 7TP~ ( n ijklm — nijk ) . r P TrP— ("ijfclm— "ijfei) , r p r rP— ("ijfci — "i jfc ) 

-"■f/l — ^DD3 U U1 > n DT3 u Ul ' lj DU3 u U\ ' ^ DT2 U Ul 

Tv- —TP jjP—{nijklm+nijk) 1 rP TjP-nij k i m jp jjP— {nijklm— nijkl) . rP jjP—{nijki+nijk) , T p jjP-riijk 
J^UB2 — lj DB3 u UB2 ' n D3 u UB2 ' ^DU3 U UB2 ' ^DB2 U UB2 ' ^D1 U UB2 

, rP TjP-^j , rP rrP-"i , TPjjP-^ 
> lj B u UB2 ' n W u UB2 "r ^S U UB2' 

TV" _ r P TT p-{nijkl m +nij k ) T p TT p-{nij k i m -nij k ) T p TT p-(nij klm -nij k i) T p TT p-nij T p jjp-rij 
JX U2 — ^ D B3 U U2 ~ r - L DT3 U U2 +lj DU3 U U2 + ^ B u U2 +lj W U U2 

1TPTTP- 1 
+- L S U U2 i 

tv- _ f P jjP-riij k i m T p TJ p-{nij k i m ~nij k ) jV TT P-(riij k im ) , r p j T P-{nijki-nijk) . T p jjP-riij 

r ^UT2 J -'D3 UT2 ' lj BT3 u UT2 ' lj DU3 u UT2 ' lj DT2 u UT2 ' lj B u UT2 

1 rP 7-rP-™ 3 , jPjjp-1 

tv- _ rP jjP-(nij k im+nij k i) , T-p jjP-(nijkl+nijk) , rP T r p ~' ni i k ' i rP 7"r p_ ("ijfei -»*«* ) i rP T T p ~ ni i k 

J^UD3 — Lj BD3 KJ VD3 ' ^DB2 U UD3 ~ rlj D2 u UD3 ' lj DT2 u UD3 ' ^D1 U UD3 

1 TP TjP- n ii 1 rP TjP- n i 1 tPttP- 1 
' t1j B u UD3 ^ ^W U UD3 ~ tlJ S u UD3> 

Lf —TP j T P-( n ijkl m +nijk) , rP TT P-(nijkl+n ijk ) jP TT p-n ijkl jP TT p-n ijk jV jjp-nij jV jjp-rij 

J^UB3 — lj DB3 u UB3 ' ^DB2 U UB3 ' ^D2 U UB3 ~< Ij Dl u UB3 ' lj B u UB3 ^ n W u UB3 

\tPtjP- 1 

tv- _ rP T jP-{nij k i+nij k ) T p TT p-{nijkl-nijk) . T p rrP-ntj T p TT p-rij T p TT p-l 
KU3 — ^BB2 U U3 +Ij DT2 U U3 +lj B U U3 +- L W U U3 +1 ^S U U3 ' 

tv- _ rP TTP-(nijkl m -ni jk ) jp jrp-n % j kl jp jjP-(riij k i -n ijk ) . T p TjP-riij jV jjP-rij T pjjp-1 
J^UT3 — n BT3 u UT3 ~ rlj D2 u UT3 ' ^ DT2 U UT3 ' lj B u UT3 ~ tIj W u UT3 ~ r2j S u UT3' 

K —TP TjP-{ n Hklm-nij k i) T p jjp-nijk , rP jjP-^ij T p jjp-nj T p TT p-l / nl \ 

KUU3 — ^DU3 U UU3 +- L D1 U UU3 +lj B U UU3 +^W U UU3 +1 ^S U UU3' \ Ui -) 
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